#include <iostream>
#include <immintrin.h>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX、AVX
#include<iomanip>
#include<vector>
#include<sys/time.h>
using namespace std;
int n ;
//float** m;
vector<float**>v;
float** generate()
{
	float** m = new float*[n];
        for(int p = 0; p < n ;p ++)
        {
                m[p] = new float[n];
                for(int q = 0; q < n; q++)
                {
                        m[p][q] = 0.0;
                }
	}
    for (int i = 0; i < n; i++)
    {
        m[i][i] = 1.0;
        for (int p = i + 1; p < n; p++)
        {
            m[i][p] = rand()%100 + 1;
        }
    }
    for (int k = 0; k < n; k++)
    {
        for (int i = k + 1; i < n; i++)
        {
            for (int p = 0; p < n; p++)
            {
                m[i][p] += m[k][p];
            }
        }
    }
    return m;
}
void init()
{
    for(int i = 0;i < 20;i++)
    {
	float** t = new float*[n];
	for(int p = 0; p < n ;p ++)
	{
		t[i] = new float[n];
		for(int q = 0; q < n; q++)
		{
			t[i][q] = 0.0;
		}
	}
	v.push_back(t);
    }
}
void show(float**m)
{
    for (int i = 0; i < n; i++)
    {
        for (int p = 0; p < n; p++)
        {
            cout << setw(5)<<m[i][p] << " ";
        }
        cout << endl;
    }
}
void gaosi(float**m)
{
    int temp = 0;
    double bei = 0;
    bool whetherFound = 0;
    int j = 0;
    __m256 t1, t2, vaki, vaij, vakj, vx;
    for (int i = 0; i < n; i++)//似乎不用考虑主对角线是0
    {
        t1 = _mm256_set1_ps(m[i][i]);
        for (j = i + 1; j + 8 < n; j += 8)
        {
            t2 = _mm256_loadu_ps(&m[i][j]);
            t2 = _mm256_div_ps(t2, t1);
            _mm256_storeu_ps(&m[i][j], t2);
        }
        for (; j < n; j++)
        {
            m[i][j] /= m[i][i];
        }
        m[i][i] = 1;
        for (int k = i + 1; k < n; k++)
        {
            vaki = _mm256_set1_ps(m[k][i]);
            for (j = i + 1; j + 8 < n; j += 8)
            {
                vaij = _mm256_loadu_ps(&m[i][j]);
                vakj = _mm256_loadu_ps(&m[k][j]);
                vx = _mm256_mul_ps(vaij, vaki);
                vakj = _mm256_sub_ps(vakj, vx);
                _mm256_storeu_ps(&m[k][j], vakj);
            }
            for (; j < n; j++)
            {
                m[k][j] -= m[i][j] * m[k][i];
            }
            m[k][i] = 0;
        }
    }
    //int temp = 0;
    //double bei = 0;
    //bool whetherFound = 0;
    //int j = 0;
    //__m512 t1, t2, vaki, vaij, vakj, vx;
    //for (int i = 0; i < n; i++)//似乎不用考虑主对角线是0
    //{
    //    t1 = _mm512_set1_ps(m[i][i]);
    //    for (j = i + 1; j + 16 < n; j += 16) 
    //    {
    //        t2 = _mm512_loadu_ps(&m[i][j]);
    //        t2 = _mm512_div_ps(t2,t1);
    //        _mm512_storeu_ps(&m[i][j], t2);
    //    }
    //    for (; j < n; j++)
    //    {
    //        m[i][j] /= m[i][i];
    //    }
    //    m[i][i] = 1;
    //    for (int k = i + 1; k < n; k++)
    //    {
    //        vaki = _mm512_set1_ps(m[k][i]);
    //        for (j = i + 1; j + 16 < n; j += 16)
    //        {
    //            vaij = _mm512_loadu_ps(&m[i][j]);
    //            vakj = _mm512_loadu_ps(&m[k][j]);
    //            vx = _mm512_mul_ps(vaij, vaki);
    //            vakj = _mm512_sub_ps(vakj, vx);
    //            _mm512_storeu_ps(&m[k][j], vakj);
    //        }
    //        for (; j < n; j++)
    //        {
    //            m[k][j] -= m[i][j] * m[k][i];
    //        }
    //        m[k][i] = 0;
    //    }
    //}
}
void commonGaosi(float**m)
{
	for(int i = 0 ;i < n; i++)
	{
        	for (int j = i + 1 ; j < n; j++)
        	{
        	    m[i][j] /= m[i][i];
        	}
		for (int k = i + 1; k < n; k++)
        	{
            	for (int j = i + 1; j < n; j++)
            	{
                	m[k][j] -= m[i][j] * m[k][i];
            	}
            	m[k][i] = 0;
        	}
	}
}
int main()
{
for(int i = 1; i <= 20; i++)
{
    n = i * 100;
    float**my = generate();
    struct timeval t1,t2;
    double timeuse;
    gettimeofday(&t1,NULL);
    gaosi(my);

    //commonGaosi(my);
    gettimeofday(&t2,NULL);
    timeuse = (t2.tv_sec - t1.tv_sec) + (double)(t2.tv_usec - t1.tv_usec)/1000000.0;
    //show(my);
    cout<<timeuse<<endl;  //输出时间（单位：ｓ）
    
}
    return 0;
}
